Write a custom CUDA kernel to optimize BPR Loss (Bayesian Personalized Ranking) with 2D inputs.

Formula: Loss = -log(sigmoid(pos - neg)) which is mathematically equivalent to log(1 + exp(neg - pos)).

Problem Analysis:
The inputs are typically 2D tensors of shape (Batch_Size, Num_Items). The operation is strictly element-wise. The standard PyTorch implementation involves chaining subtraction, logsigmoid (which implies exp, add, log, reciprocal), and potentially a final reduction. This chain creates intermediate tensors and consumes significant memory bandwidth relative to the low arithmetic intensity.

Optimization Strategy: Fused Element-wise Kernel with Vectorization

1. Flattened Processing: Since the operation is element-wise, the kernel treats the 2D contiguous input tensors as flattened 1D arrays. This simplifies indexing logic and allows for uniform block distribution regardless of specific 2D dimensions.

2. Vectorized Loads (float4): The kernel uses float4 data types to load and store 4 float elements (128 bits) per instruction. This drastically reduces the number of memory transactions, which is critical for memory-bound operations like this.

3. Fused In-Register Math:
   - Load 4 elements from pos_scores and neg_scores.
   - Compute difference: x = neg - pos.
   - Apply stable Softplus function: if x > 0 return x + log(1 + exp(-x)), else return log(1 + exp(x)).
   - Store the result.

4. Reduction Handling: The kernel calculates element-wise losses. The final reduction (mean or sum) is handled by the C++ wrapper using optimized ATen primitives, avoiding the overhead of returning to Python for aggregation.
  
Here's an example to show you the syntax of inline embedding custom CUDA operators in torch: The example given architecture is:   
  
```python
import torch
import torch.nn as nn
import torch.nn.functional as F

BATCH_SIZE = 16384
NUM_ITEMS = 1024
SHAPE = (BATCH_SIZE, NUM_ITEMS)

class BPRLoss(nn.Module):
    def __init__(self, reduction='mean'):
        super(BPRLoss, self).__init__()
        self.reduction = reduction

    def forward(self, pos_scores: torch.Tensor, neg_scores: torch.Tensor) -> torch.Tensor:
        diff = pos_scores - neg_scores

        loss = -F.logsigmoid(diff)
        
        if self.reduction == 'mean':
            return loss.mean()
        elif self.reduction == 'sum':
            return loss.sum()
        return loss

class Model(nn.Module):
    def __init__(self, reduction='none'):
        super(Model, self).__init__()
        self.loss_fn = BPRLoss(reduction=reduction)
    
    def forward(self, pos, neg):
        return self.loss_fn(pos, neg)

def get_inputs():
    pos_scores = torch.randn(SHAPE, dtype=torch.float32)
    neg_scores = torch.randn(SHAPE, dtype=torch.float32)
    return [pos_scores.contiguous(), neg_scores.contiguous()]

def get_init_inputs():
    return ['none']